tmap

Dependências do sistema operacional para a biblioteca tmap:

sudo apt install libudunits2-dev libgdal-dev libgeos-dev libproj-dev libfontconfig1-dev

base de dados que criaremos um ranking com análise fatorial:

library(textshape)
data_original = read.csv("data/atlas.csv")
data_original = column_to_rownames(data_original,'distritos')
data = data_original[,-1]

Correlações das colunas numéricas:

library(PerformanceAnalytics)
chart.Correlation(data)

Teste de esfericidade de bartlet

library(psych)
rho = cor(data)
cortest.bartlett(R = rho)

## $chisq
## [1] 780.9853
## 
## $p.value
## [1] 8.646135e-141
## 
## $df
## [1] 36

Como p-value é 8.646135e-141 dizemos que na significância de 5% (ou nível confiança de 95%) ficamos com a hipótese alternativa de que essa matriz de correlação é estatisticamente diferente da matriz de identidade (zerada)

Normanlizando base de dados:

data_std = data.frame(scale(data))

Modelo:

modelo = prcomp(data_std)
summary(modelo)

## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.2262 1.0790 0.9982 0.85092 0.72753 0.63113 0.36010
## Proportion of Variance 0.5507 0.1294 0.1107 0.08045 0.05881 0.04426 0.01441
## Cumulative Proportion  0.5507 0.6800 0.7907 0.87120 0.93001 0.97427 0.98868
##                            PC8     PC9
## Standard deviation     0.25508 0.19196
## Proportion of Variance 0.00723 0.00409
## Cumulative Proportion  0.99591 1.00000

Autovalores:

modelo$sdev^2

## [1] 4.95603069 1.16433814 0.99635412 0.72406663 0.52930223 0.39832415 0.12966871
## [8] 0.06506649 0.03684884

O primeiro fator captura variância de quase 5 variáves.

Autovetores:

autovetores = modelo$rotation

Autovetor, ou seja, o peso da cada variável em cada fator:

library(reshape2)
data_melt=melt(autovetores)
colnames(data_melt) = c('Variavel','Fator','value')

library(tidyverse)
ggplot(aes(x = Variavel, y = value, fill = Variavel),data=data_melt) +
  geom_bar(stat = "identity") + 
  facet_wrap(~Fator) +
  theme(axis.text.x = element_text(angle = 90))

Com o screeplot podemos ver o decaimento da variância compartilhada:

library(factoextra)
fviz_eig(modelo,addlabels=T)

cargas fatorais (correlação entre os fatores estabelecidos e as variáveis originais):

cargas_fatoriais = modelo$rotation %*% diag(modelo$sdev)
cargas_fatoriais

##                    [,1]        [,2]         [,3]        [,4]         [,5]
## renda        -0.8321269  0.36897098 -0.172642067  0.05394096 -0.248818640
## quota        -0.9006058  0.22645316 -0.154731604  0.11909773 -0.202225206
## escolaridade -0.9665432 -0.02459367 -0.009340238  0.07197491 -0.002825683
## idade        -0.9601351 -0.06544713  0.063367386  0.10798882  0.022849354
## mortalidade   0.6556993 -0.17663396  0.173477851  0.49370067 -0.495545240
## txcresc       0.6967815  0.33744591  0.277249011 -0.38556544 -0.254712046
## causasext     0.6665775 -0.04390999 -0.438317349  0.41725042  0.210756174
## favel         0.4571854  0.44283981 -0.686884484 -0.13511980 -0.106507880
## denspop      -0.1662813 -0.79304783 -0.409402407 -0.32461428 -0.244460635
##                     [,6]        [,7]         [,8]          [,9]
## renda         0.19387656 -0.12866888  0.144109283 -0.0434546573
## quota         0.13422152 -0.05322146 -0.188989739  0.0451875103
## escolaridade -0.03266914  0.19137260  0.084723201  0.1226578170
## idade        -0.04106977  0.19269177 -0.032852846 -0.1332741152
## mortalidade  -0.13537427  0.03179396  0.011028723  0.0004999632
## txcresc       0.29475000  0.15232212 -0.011818498 -0.0007991541
## causasext     0.37010628  0.07805418  0.005643949 -0.0022772604
## favel        -0.29519474  0.07887499 -0.004964341 -0.0087769279
## denspop       0.10317107  0.00055644  0.002596192 -0.0053459708

Mantendo somente os dois primeiros fatores:

cargas_fatoriais = cargas_fatoriais[,1:2]
cargas_fatoriais

##                    [,1]        [,2]
## renda        -0.8321269  0.36897098
## quota        -0.9006058  0.22645316
## escolaridade -0.9665432 -0.02459367
## idade        -0.9601351 -0.06544713
## mortalidade   0.6556993 -0.17663396
## txcresc       0.6967815  0.33744591
## causasext     0.6665775 -0.04390999
## favel         0.4571854  0.44283981
## denspop      -0.1662813 -0.79304783

Comunalidade, variâncias totais compartilhadas em cada um dos fatores considerados:

comunalidade = rowSums(cargas_fatoriais^2)
comunalidade

##        renda        quota escolaridade        idade  mortalidade      txcresc 
##    0.8285748    0.8623719    0.9348107    0.9261428    0.4611411    0.5993741 
##    causasext        favel      denspop 
##    0.4462537    0.4051255    0.6565743

Capturamos 82,86% da variância de renda e apenas 40,5% de favel, considerando apenas dois fatores.

é importante mostrar as cargas_fatoriais junto com a comunalidade:

cbind(cargas_fatoriais,comunalidade)

##                                     comunalidade
## renda        -0.8321269  0.36897098    0.8285748
## quota        -0.9006058  0.22645316    0.8623719
## escolaridade -0.9665432 -0.02459367    0.9348107
## idade        -0.9601351 -0.06544713    0.9261428
## mortalidade   0.6556993 -0.17663396    0.4611411
## txcresc       0.6967815  0.33744591    0.5993741
## causasext     0.6665775 -0.04390999    0.4462537
## favel         0.4571854  0.44283981    0.4051255
## denspop      -0.1662813 -0.79304783    0.6565743

Gráfico das cargas fatoriais:

library(tidyverse)
library(ggrepel)

round(summary(modelo)$importance[2,1] * 100,digits = 2)

## [1] 55.07

round(summary(modelo)$importance[2,2] * 100,digits = 2)

## [1] 12.94

cargas_fatoriais = data.frame(cargas_fatoriais)

ggplot(data=cargas_fatoriais,
  aes(x = X1, y = X2)) +
  geom_point(color = "dodgerblue4") +
  geom_hline(yintercept = 0, color = "darkgoldenrod3", linetype = "dashed") +
  geom_vline(xintercept = 0, color = "darkgoldenrod3", linetype = "dashed") +
  geom_text_repel(label = row.names(cargas_fatoriais)) +
  labs(x = 'F1 (62,98%)',y = 'F2 (25,01%)')

densipop se relaciona mais com F2 e as demais se relacionam todas com F1.

scores_fatoriais:

scores_fatoriais = t(modelo$rotation)/modelo$sdev
scores_fatoriais = t(scores_fatoriais)
scores_fatoriais

##                      PC1         PC2          PC3         PC4          PC5
## renda        -0.16790189  0.31689332 -0.173273803  0.07449723 -0.470088022
## quota        -0.18171917  0.19449089 -0.155297801  0.16448449 -0.382059990
## escolaridade -0.19502366 -0.02112244 -0.009374416  0.09940371 -0.005338506
## idade        -0.19373066 -0.05620973  0.063599261  0.14914211  0.043168823
## mortalidade   0.13230331 -0.15170332  0.174112644  0.68184425 -0.936223595
## txcresc       0.14059264  0.28981780  0.278263526 -0.53249994 -0.481222313
## causasext     0.13449826 -0.03771240 -0.439921248  0.57625970  0.398177376
## favel         0.09224829  0.38033608 -0.689397945 -0.18661238 -0.201223182
## denspop      -0.03355130 -0.68111471 -0.410900501 -0.44832100 -0.461854530
##                      PC6          PC7         PC8         PC9
## renda         0.48673063 -0.992289382  2.21480047 -1.17926809
## quota         0.33696555 -0.410441804 -2.90456349  1.22629409
## escolaridade -0.08201647  1.475858028  1.30210199  3.32867543
## idade        -0.10310639  1.486031391 -0.50491194 -3.61677946
## mortalidade  -0.33985956  0.245193798  0.16949929  0.01356795
## txcresc       0.73997522  1.174702242 -0.18163726 -0.02168736
## causasext     0.92915851  0.601950805  0.08674126 -0.06180006
## favel        -0.74109174  0.608280839 -0.07629643 -0.23818738
## denspop       0.25901285  0.004291244  0.03990061 -0.14507841

mantendo 2 scores fatoriais:

scores_fatoriais = scores_fatoriais[,1:2]
scores_fatoriais

##                      PC1         PC2
## renda        -0.16790189  0.31689332
## quota        -0.18171917  0.19449089
## escolaridade -0.19502366 -0.02112244
## idade        -0.19373066 -0.05620973
## mortalidade   0.13230331 -0.15170332
## txcresc       0.14059264  0.28981780
## causasext     0.13449826 -0.03771240
## favel         0.09224829  0.38033608
## denspop      -0.03355130 -0.68111471

Para criação do indicador manteremos apenas o primeiro fator.

Multiplicação entre primeiro score fatorial e o valores observado padronizado:

multiplicacoes_F1 <- t(apply(data_std, 1, function(x) x * scores_fatoriais[,1]))
F1 = data.frame(-rowSums(multiplicacoes_F1)) # deixamos negativo para o ranking ficar "bunitinho" - quanto maior, signica melhor, e não o contrario

Colocando fator 1 no dataframe original:

data_original["fator1"] <- F1

Criando uma coluna pontuação do ranking pela soma ponderada dos fatores por sua variância compartilhada, ou seja, vamos considerar 55,06% que é o que esse valor explica:

var_compartilhada <- (modelo$sdev ^ 2/sum(modelo$sdev ^ 2))
data_original$pontuacao = data_original$fator1 * var_compartilhada[1] # 55,06%

Abrindo o shapefile da cidade de são Paulo:

library(maptools)
sao_paulo = readShapeSpatial("data/shapefile_saopaulo/sao_paulo")

Salvando shapefile se necessário:

library(maptools)
writeSpatialShape(sao_paulo, "data/shapefile_saopaulo/sao_paulo")

Plotando mapa:

library(tmap)
tm_shape(sao_paulo) + tm_borders()

Colocando dataframe com Fator1 no mapa:

sao_paulo@data$COD_DIST = as.numeric(sao_paulo@data$COD_DIST)
data_original$distritos = rownames(data_original)

distritos_dados <- merge(sao_paulo,
                         data_original,
                         by.x = "COD_DIST",
                         by.y = "cod_ibge")

Vendo o mapa:

#tmap_mode("plot")
tmap_mode("view")

tm_shape(distritos_dados) +
  tm_fill("pontuacao", midpoint = 0, palette = "RdBu", 
          style = "quantile", n = 10, legend.show = T) +
  tm_borders(alpha = 0.8) +
  tm_text("distritos")